##################################################################################################################
# Code for Main Results 
##################################################################################################################
write.table(data.frame(Description=c("","","","Main Results")),row.names=F,col.names=F,quote=F)

##################################################################################################################
# To run, please set the working directory to match the location of the folder containing all replication files.
##################################################################################################################
#setwd("Replication_PSRM_ST24")

#if not yet created, create "output" folder
if("output"%in%list.files()==F){dir.create("output")}

#Information about R Version used
as.data.frame(R.Version())
#   platform                arch       os           system status major minor year month day svn.rev language               version.string  nickname
# 1 x86_64-apple-darwin20 x86_64 darwin20 x86_64, darwin20            4   4.0 2024    04  24   86474        R R version 4.4.0 (2024-04-24) Puppy Cup

#Clear environment
#rm(list=ls());gc();gc();gc();gc()

#If not yet installed, install package
if(length(find.package("dplyr",quiet=T))==0){install.packages("dplyr")}

###packages
library(dplyr)
#packageVersion("dplyr") #[1] '1.1.4'

#Set POSIX locale
Sys.setlocale(locale="C")

##################
#Load main dataset
##################
loball5ciA<-readRDS("Data/PSRM_main_dataset.rds")

write.table(data.frame(Description=c("","Description in Data Section:")),row.names=F,col.names=F,quote=F)
write.table(data.frame(Description=c("","Number of observations of lobbyist-client relationships:")),row.names=F,col.names=F,quote=F)
print(nrow(loball5ciA)-length(which(loball5ciA$type=="nonlobbyist"))) #43817
write.table(data.frame(Description=c("","Breakdown of observations, including non-licensed lobbyists (non-lobbyists):")),row.names=F,col.names=F,quote=F)
print(table(loball5ciA$type))
# contract     inhouse nonlobbyist 
# 28541       15276        1956 
####################
# Prepare Analyses
####################

#####
#within dyad deviations

#contract and in-house
agg67c<-loball5ciA %>%
  filter(type%in%c("contract","inhouse")==T) %>%
  group_by(principal_cleaned,lobname_cleaned, halfyear,type,session,both_ever_count1,
           inhouse_session_count1,contract_session_count1,organization_type) %>% 
  summarise(.groups="keep",sum_hours=sum(hours,na.rm=T),sum_amount=sum(Rcomparable_total_amounts_cleaned_exp,na.rm=T))

agg67d2<-agg67c
#create hours worked variable so that later log transformation puts zero weight on extensive margin
agg67d2$sum_hours_cr0<-agg67d2$sum_hours
agg67d2$sum_hours_cr0[agg67d2$sum_hours_cr0>0]<-agg67d2$sum_hours_cr0[agg67d2$sum_hours_cr0>0]/min(agg67d2$sum_hours_cr0[agg67d2$sum_hours_cr0>0])
agg67d2$sum_hours_cr0[agg67d2$sum_hours_cr0==0]<-1

length(which(agg67d2$sum_hours==0&agg67d2$sum_amount==0)) #9207

write.table(data.frame(Description=c("","Description in Empirical Approach Section:")),row.names=F,col.names=F,quote=F)
write.table(data.frame(Description=c("","Footnote 16 (and Appendix A4): Percentage of observations with zero expenditures that have non-zero hours associated with them:")),row.names=F,col.names=F,quote=F)
print(round(length(which(agg67d2$sum_hours==0&agg67d2$sum_amount==0))/length(which(agg67d2$sum_amount==0)),2)) #90%

#count number of observations
n671<-nrow(agg67d2) #[1] 43051
table(agg67d2$type) 
# contract  inhouse 
# 27775    15276

#create variables for analyses
agg67d2 <-agg67d2 %>%  
  filter(sum_amount>0) %>% 
  group_by(principal_cleaned,lobname_cleaned,type,session) %>%
  mutate(
         dif_log_mean_amount=log(sum_amount)-log(mean(sum_amount,na.rm=T)),
         dif_log_mean_hours_cr0=log(sum_hours_cr0)-log(mean(sum_hours_cr0,na.rm=T))
         )

n672<-nrow(agg67d2) #[1] 32813
#43051-32813=10238
#10238/43051=0.24
table(agg67d2$type) 
# contract  in-house 
# 20254    12559

agg67d2$count<-1
agg67d2 <-agg67d2 %>% 
  group_by(principal_cleaned,lobname_cleaned,type,session) %>%
  mutate(num_dyad_obs=sum(count)) 
#exclude observations when num_dyad_obs==1
agg67d2<-agg67d2[which(agg67d2$num_dyad_obs>1),]
#count obs
n673<-nrow(agg67d2) #[1] 30888
#32813-30888=1925
#1925/43051=0.04
table(agg67d2$type) 
# contract  inhouse 
# 18890    11998

#12163/43051=28% total removed (10238+1925)

###contract, half-yearly
agg6c<-loball5ciA %>%
  filter(type%in%c("contract")==T) %>%
  group_by(principal_cleaned,lobname_cleaned, halfyear,session,both_ever_count1,
           inhouse_session_count1,contract_session_count1,organization_type) %>% 
  summarise(.groups="keep",sum_hours=sum(hours,na.rm=T),sum_amount=sum(Rcomparable_total_amounts_cleaned_exp,na.rm=T))

#get total amount of compensation (by dyad) in biennium
agg6c <- agg6c %>% 
  group_by(lobname_cleaned,halfyear) %>%
  mutate(count_client=length(unique(principal_cleaned)))

agg6d2<-agg6c
#create hours worked variable so that later log transformation puts zero weight on extensive margin
agg6d2$sum_hours_cr0<-agg6d2$sum_hours
agg6d2$sum_hours_cr0[agg6d2$sum_hours_cr0>0]<-agg6d2$sum_hours_cr0[agg6d2$sum_hours_cr0>0]/min(agg6d2$sum_hours_cr0[agg6d2$sum_hours_cr0>0])
agg6d2$sum_hours_cr0[agg6d2$sum_hours_cr0==0]<-1

n61<-nrow(agg6d2) # 27775

agg6d2 <-agg6d2 %>%  
  filter(sum_amount>0) %>% 
  group_by(principal_cleaned,lobname_cleaned,session) %>%
  mutate(
         dif_log_mean_amount=log(sum_amount)-log(mean(sum_amount,na.rm=T)),
         dif_log_mean_hours_cr0=log(sum_hours_cr0)-log(mean(sum_hours_cr0,na.rm=T))
         )

#count #obs
n62<-nrow(agg6d2) #20254
#27775-20254=7521
#7521/27775=0.27
agg6d2$count<-1
agg6d2 <-agg6d2 %>% 
  group_by(principal_cleaned,lobname_cleaned,session) %>%
  mutate(num_dyad_obs=sum(count)) 
#exclude observations when num_dyad_obs==1
agg6d2<-agg6d2[which(agg6d2$num_dyad_obs>1),]
#count obs
n63<-nrow(agg6d2) #18890
#20254-18890= 1364
# 1364/27775=0.05

#####in-house

#in-house
agg7c<-loball5ciA %>%
  filter(type%in%c("inhouse")==T)%>%
  group_by(principal_cleaned,lobname_cleaned,halfyear,session,both_ever_count1,
           inhouse_session_count1,contract_session_count1,organization_type) %>% 
  summarise(.groups="keep",sum_hours=sum(hours,na.rm=T),sum_amount=sum(Rcomparable_total_amounts_cleaned_exp,na.rm=T))

#get number of clients for whom lobbyist is working in halfyear
agg7c <- agg7c %>% group_by(lobname_cleaned,halfyear) %>%
  mutate(count_client=length(unique(principal_cleaned)))



agg7d2<-agg7c
#create hours worked variable so that later log transformation puts zero weight on extensive margin
agg7d2$sum_hours_cr0<-agg7d2$sum_hours
agg7d2$sum_hours_cr0[agg7d2$sum_hours_cr0>0]<-agg7d2$sum_hours_cr0[agg7d2$sum_hours_cr0>0]/min(agg7d2$sum_hours_cr0[agg7d2$sum_hours_cr0>0])
agg7d2$sum_hours_cr0[agg7d2$sum_hours_cr0==0]<-1

n71<-nrow(agg7d2) #15276

#create variables for analyses
agg7d2 <-agg7d2 %>%  
  filter(sum_amount>0) %>% 
  group_by(principal_cleaned,lobname_cleaned,session) %>%
  mutate(
         dif_log_mean_amount=log(sum_amount)-log(mean(sum_amount,na.rm=T)),
         dif_log_mean_hours_cr0=log(sum_hours_cr0)-log(mean(sum_hours_cr0,na.rm=T))
         )

n72<-nrow(agg7d2) # 12559
# 15276-12559= 2717
# 2717/15276=0.18

agg7d2$count<-1
agg7d2 <-agg7d2 %>% group_by(principal_cleaned,lobname_cleaned,session) %>%
  mutate(num_dyad_obs=sum(count)) 
#exclude observations when num_dyad_obs==1
agg7d2<-agg7d2[which(agg7d2$num_dyad_obs>1),]

n73<-nrow(agg7d2) # 11998
# 12559-11998=561
# 561/15276=0.04

###examine prevalence of observations with zero hours (footnote 17)
write.table(data.frame(Description=c("","Footnote 17: Percentage of observations in the analysis sample that have zero hours worked:")),row.names=F,col.names=F,quote=F)
print(round(length(which(agg67d2$sum_hours==0))/nrow(agg67d2),3)) #0.004
length(which(agg6d2$sum_hours==0))/nrow(agg6d2) #0.0046
length(which(agg7d2$sum_hours==0))/nrow(agg7d2) #0.0019

sum(agg67d2$sum_amount[which(agg67d2$sum_hours==0)])/sum(agg67d2$sum_amount) #0.0039
sum(agg6d2$sum_amount[which(agg6d2$sum_hours==0)])/sum(agg6d2$sum_amount) #0.0060
sum(agg7d2$sum_amount[which(agg7d2$sum_hours==0)])/sum(agg7d2$sum_amount) #0.0008


#see footnote 18: (3% of the in-house lobbyist observations in the analysis sample are from lobbyists that work for multiple clients)
write.table(data.frame(Description=c("","Footnote 18: Percentage of the in-house lobbyist observations in the analysis sample are from lobbyists that work for multiple clients:")),row.names=F,col.names=F,quote=F)
print(round(length(which(agg7c$count_client>1))/nrow(agg7c),3)) #0.03 

#see footnote 18: 4.3% of contract lobbyists observations are from lobbyists that work for a single client at a given time
write.table(data.frame(Description=c("","Footnote 18: Percentage of contract lobbyists observations are from lobbyists that work for a single client at a given time:")),row.names=F,col.names=F,quote=F)
print(round(length(which(agg6c$count_client==1))/nrow(agg6c),3)) #0.043



############################################################
############ Aggregated to yearly
############################################################

#contract and in-house, yearly
agg67d3<-agg67c %>%
  group_by(principal_cleaned,lobname_cleaned,year=substr(halfyear,1,4),session, type,both_ever_count1,
           inhouse_session_count1,contract_session_count1,organization_type) %>%
  summarise(.groups="keep",sum_hours=sum(sum_hours,na.rm=T),sum_amount=sum(sum_amount,na.rm=T))

#create hours worked variable so that later log transformation puts zero weight on extensive margin
agg67d3$sum_hours_cr0<-agg67d3$sum_hours
agg67d3$sum_hours_cr0[agg67d3$sum_hours_cr0>0]<-agg67d3$sum_hours_cr0[agg67d3$sum_hours_cr0>0]/min(agg67d3$sum_hours_cr0[agg67d3$sum_hours_cr0>0])
agg67d3$sum_hours_cr0[agg67d3$sum_hours_cr0==0]<-1

n67a1<-nrow(agg67d3) #22611

#create variables for analyses
agg67d3 <-agg67d3 %>%  
  filter(sum_amount>0) %>% 
  group_by(principal_cleaned,lobname_cleaned,type) %>%
  mutate(
         dif_log_mean_amount=log(sum_amount)-log(mean(sum_amount,na.rm=T)),
         dif_log_mean_hours_cr0=log(sum_hours_cr0)-log(mean(sum_hours_cr0,na.rm=T))
         )

n67a2<-nrow(agg67d3) #18902
#22611-18902=3709
#3709/22611=0.16
agg67d3$count<-1

agg67d3 <-agg67d3 %>% 
  group_by(principal_cleaned,lobname_cleaned,type) %>%
  mutate(num_dyad_obs=sum(count)) 
#exclude observations when num_dyad_obs==1
agg67d3<-agg67d3[which(agg67d3$num_dyad_obs>1),]

n67a3<-nrow(agg67d3) #17404
#18902-17404=1498
#1498/22611=0.07

agg67d3<-agg67d3[order(agg67d3$principal_cleaned,agg67d3$lobname_cleaned,as.numeric(agg67d3$year),agg67d3$type),]
agg67d3<-agg67d3 %>% group_by(principal_cleaned,lobname_cleaned,type) %>%
  mutate(nyears=cumsum(count))

###contract lobbyists
agg6d3<-agg6c %>%
  group_by(principal_cleaned,lobname_cleaned,year=substr(halfyear,1,4),session, both_ever_count1,
           inhouse_session_count1,contract_session_count1,organization_type) %>%
  summarise(.groups="keep",sum_hours=sum(sum_hours,na.rm=T),sum_amount=sum(sum_amount,na.rm=T))

n6a1<-nrow(agg6d3) #14594

#create hours worked variable so that later log transformation puts zero weight on extensive margin
agg6d3$sum_hours_cr0<-agg6d3$sum_hours
agg6d3$sum_hours_cr0[agg6d3$sum_hours_cr0>0]<-agg6d3$sum_hours_cr0[agg6d3$sum_hours_cr0>0]/min(agg6d3$sum_hours_cr0[agg6d3$sum_hours_cr0>0])
agg6d3$sum_hours_cr0[agg6d3$sum_hours_cr0==0]<-1

#create variables for analyses
agg6d3 <-agg6d3 %>%  
  filter(sum_amount>0) %>% 
  group_by(principal_cleaned,lobname_cleaned) %>%
  mutate(
    dif_log_mean_amount=log(sum_amount)-log(mean(sum_amount,na.rm=T)),
    dif_log_mean_hours_cr0=log(sum_hours_cr0)-log(mean(sum_hours_cr0,na.rm=T))
    )

n6a2<-nrow(agg6d3) # 11876
#14594-11876=2718
#2718/14594=0.11
agg6d3$count<-1

agg6d3 <-agg6d3 %>% group_by(principal_cleaned,lobname_cleaned) %>%
  mutate(num_dyad_obs=sum(count)) 

#exclude observations when num_dyad_obs==1
agg6d3<-agg6d3[which(agg6d3$num_dyad_obs>1),]
n6a3<-nrow(agg6d3) # 10751
#11876-10751=1125
#1125/14594=0.08

agg6d3<-agg6d3[order(agg6d3$principal_cleaned,agg6d3$lobname_cleaned,as.numeric(agg6d3$year)),]
agg6d3<-agg6d3 %>% group_by(principal_cleaned,lobname_cleaned) %>%
  mutate(nyears=cumsum(count))

#
agg7d3<-agg7c %>%
  group_by(principal_cleaned,lobname_cleaned,year=substr(halfyear,1,4),session, both_ever_count1,
           inhouse_session_count1,contract_session_count1,organization_type) %>%
  summarise(.groups="keep",sum_hours=sum(sum_hours,na.rm=T),sum_amount=sum(sum_amount,na.rm=T))

n7a1<-nrow(agg7d3) # 8017

#create hours worked variable so that later log transformation puts zero weight on extensive margin
agg7d3$sum_hours_cr0<-agg7d3$sum_hours
agg7d3$sum_hours_cr0[agg7d3$sum_hours_cr0>0]<-agg7d3$sum_hours_cr0[agg7d3$sum_hours_cr0>0]/min(agg7d3$sum_hours_cr0[agg7d3$sum_hours_cr0>0])
agg7d3$sum_hours_cr0[agg7d3$sum_hours_cr0==0]<-1

#create variables for analyses
agg7d3 <-agg7d3 %>%  
  filter(sum_amount>0) %>% 
  group_by(principal_cleaned,lobname_cleaned) %>%
  mutate(
    dif_log_mean_amount=log(sum_amount)-log(mean(sum_amount,na.rm=T)),
    dif_log_mean_hours_cr0=log(sum_hours_cr0)-log(mean(sum_hours_cr0,na.rm=T))
    )

n7a2<-nrow(agg7d3) #7026
#8017-7026=991
#991/8017=0.06
agg7d3$count<-1

agg7d3 <-agg7d3 %>% group_by(principal_cleaned,lobname_cleaned) %>%
  mutate(num_dyad_obs=sum(count)) 

#exclude observations when num_dyad_obs==1
agg7d3<-agg7d3[which(agg7d3$num_dyad_obs>1),]
n7a3<-nrow(agg7d3) #6653
#7026-6653=373
#373/8017 5%

agg7d3<-agg7d3[order(agg7d3$principal_cleaned,agg7d3$lobname_cleaned,as.numeric(agg7d3$year)),]
agg7d3<-agg7d3 %>% group_by(principal_cleaned,lobname_cleaned) %>%
  mutate(nyears=cumsum(count))

png("output/PSRM_fig1.png",width=12,height=8.83,units="in", res=300)
  par(mfrow=c(2,3))
  #Panel 1
  tt<-cor.test(agg67d2$dif_log_mean_amount,agg67d2$dif_log_mean_hours_cr0)
  plot(main=c("(1) All Lobbyists"),cex.main=1.5,cex.lab=1.2,cex.axis=1.5, x=agg67d2$dif_log_mean_amount,y=agg67d2$dif_log_mean_hours_cr0,
       xlim=c(-8,2),ylim=c(-8,2),cex=1,col="gray",xlab="Ln(Dyad Amount)-Ln(Dyad Biennium Mean Amount)",ylab="Ln(Dyad Hours)-Ln(Dyad Biennium Mean Hours)")
  text(x=-5.75,y=1.85,paste("r = ",round(tt$estimate,3)," (",round(sqrt( (1-tt$estimate^2)/(tt$parameter)),3),")",sep=""),cex=1.5)
  legend(bg="white",bty = "n","bottomright",paste("n = ",floor(tt$parameter/1000),",", 2+tt$parameter%%1000,sep=""),cex=1.5)
  lines(lowess(x=agg67d2$dif_log_mean_amount,y=agg67d2$dif_log_mean_hours_cr0,f=0.68),lty=2)
  #Panel 2
  tt<-cor.test(agg6d2$dif_log_mean_amount,agg6d2$dif_log_mean_hours_cr0)
  plot(main=c("(2) Contract Lobbyists"),cex.main=1.5,cex.lab=1.2,cex.axis=1.5,x=agg6d2$dif_log_mean_amount,y=agg6d2$dif_log_mean_hours_cr0,
       xlim=c(-8,2),ylim=c(-8,2),cex=1,col="gray",xlab="Ln(Dyad Amount)-Ln(Dyad Biennium Mean Amount)",ylab="Ln(Dyad Hours)-Ln(Dyad Biennium Mean Hours)")
  text(x=-5.75,y=1.85,paste("r = ",round(tt$estimate,3)," (",round(sqrt( (1-tt$estimate^2)/(tt$parameter)),3),")",sep=""),cex=1.5)
  legend(bg="white",bty = "n","bottomright",paste("n = ",floor(tt$parameter/1000),",", 2+tt$parameter%%1000,sep=""),cex=1.5)
  lines(lowess(x=agg6d2$dif_log_mean_amount,y=agg6d2$dif_log_mean_hours_cr0,f=0.68) ,lty=2)
  #Panel 3
  tt<-cor.test(agg7d2$dif_log_mean_amount,agg7d2$dif_log_mean_hours_cr0)
  plot(main=c("(3) In-House Lobbyists"),cex.main=1.5,cex.lab=1.2,cex.axis=1.5,x=agg7d2$dif_log_mean_amount,y=agg7d2$dif_log_mean_hours_cr0,
       xlim=c(-8,2),ylim=c(-8,2),cex=1,col="gray",xlab="Ln(Dyad Amount)-Ln(Dyad Biennium Mean Amount)",ylab="Ln(Dyad Hours)-Ln(Dyad Biennium Mean Hours)")
  text(x=-5.75,y=1.85,paste("r = ",round(tt$estimate,3)," (",round(sqrt( (1-tt$estimate^2)/(tt$parameter)),3),")",sep=""),cex=1.5)
  legend(bg="white",bty = "n","bottomright",paste("n = ",floor(tt$parameter/1000),",", 2+tt$parameter%%1000,sep=""),cex=1.5)
  lines(lowess(x=agg7d2$dif_log_mean_amount,y=agg7d2$dif_log_mean_hours_cr0,f=0.68) ,lty=2)
  #Panel 4
  tt<-cor.test(agg67d3$dif_log_mean_amount,agg67d3$dif_log_mean_hours_cr0)
  plot(main=c("(4) All Lobbyists - Yearly Aggregation"),cex.main=1.5,cex.lab=1.2,cex.axis=1.5,x=agg67d3$dif_log_mean_amount,y=agg67d3$dif_log_mean_hours_cr0,
       xlim=c(-8,2),ylim=c(-8,2),cex=1,col="gray",xlab="Ln(Dyad Amount)-Ln(Dyad Mean Amount)",ylab="Ln(Dyad Hours)-Ln(Dyad Mean Hours)")
  text(x=-5.75,y=1.85,paste("r = ",round(tt$estimate,3)," (",round(sqrt( (1-tt$estimate^2)/(tt$parameter)),3),")",sep=""),cex=1.5)
  legend(bg="white",bty = "n","bottomright",paste("n = ",floor(tt$parameter/1000),",", 2+tt$parameter%%1000,sep=""),cex=1.5)
  lines(lowess(x=agg67d3$dif_log_mean_amount,y=agg67d3$dif_log_mean_hours_cr0,f=0.68) ,lty=2)
  #Panel 5
  tt<-cor.test(agg6d3$dif_log_mean_amount,agg6d3$dif_log_mean_hours_cr0)
  plot(main=c("(5) Contract Lobbyists - Yearly Aggregation"),cex.main=1.5,cex.lab=1.2,cex.axis=1.5,x=agg6d3$dif_log_mean_amount,y=agg6d3$dif_log_mean_hours_cr0,
       xlim=c(-8,2),ylim=c(-8,2),cex=1,col="gray",xlab="Ln(Dyad Amount)-Ln(Dyad Mean Amount)",ylab="Ln(Dyad Hours)-Ln(Dyad Mean Hours)")
  text(x=-5.75,y=1.85,paste("r = ",round(tt$estimate,3)," (",round(sqrt( (1-tt$estimate^2)/(tt$parameter)),3),")",sep=""),cex=1.5)
  legend(bg="white",bty = "n","bottomright",paste("n = ",floor(tt$parameter/1000),",", 2+tt$parameter%%1000,sep=""),cex=1.5)
  lines(lowess(x=agg6d3$dif_log_mean_amount,y=agg6d3$dif_log_mean_hours_cr0,f=0.68) ,lty=2)
  #Panel 6
  tt<-cor.test(agg7d3$dif_log_mean_amount,agg7d3$dif_log_mean_hours_cr0)
  plot(main=c("(6) In-House Lobbyists - Yearly Aggregation"),cex.main=1.5,cex.lab=1.2,cex.axis=1.5,x=agg7d3$dif_log_mean_amount,y=agg7d3$dif_log_mean_hours_cr0,
       xlim=c(-8,2),ylim=c(-8,2),cex=1,col="gray",xlab="Ln(Dyad Amount)-Ln(Dyad Mean Amount)",ylab="Ln(Dyad Hours)-Ln(Dyad Mean Hours)")
  text(x=-5.75,y=1.85,paste("r = ",round(tt$estimate,3)," (",round(sqrt( (1-tt$estimate^2)/(tt$parameter)),3),")",sep=""),cex=1.5)
  legend(bg="white",bty = "n","bottomright",paste("n = ",floor(tt$parameter/1000),",", 2+tt$parameter%%1000,sep=""),cex=1.5)
  lines(lowess(x=agg7d3$dif_log_mean_amount,y=agg7d3$dif_log_mean_hours_cr0,f=0.68) ,lty=2)

dev.off()

write.table(data.frame(Description=c("(N.B. The file for this figure is saved as 'PSRM_fig1.png' in the 'output' folder.)")),row.names=F,col.names=F,quote=F)

#Include code again so that Figure A6 also shown in R if code file is run using source()
par(mfrow=c(2,3))
#Panel 1
tt<-cor.test(agg67d2$dif_log_mean_amount,agg67d2$dif_log_mean_hours_cr0)
plot(main=c("(1) All Lobbyists"),cex.main=1.5,cex.lab=1.2,cex.axis=1.5, x=agg67d2$dif_log_mean_amount,y=agg67d2$dif_log_mean_hours_cr0,
     xlim=c(-8,2),ylim=c(-8,2),cex=1,col="gray",xlab="Ln(Dyad Amount)-Ln(Dyad Biennium Mean Amount)",ylab="Ln(Dyad Hours)-Ln(Dyad Biennium Mean Hours)")
text(x=-5.75,y=1.85,paste("r = ",round(tt$estimate,3)," (",round(sqrt( (1-tt$estimate^2)/(tt$parameter)),3),")",sep=""),cex=1.5)
legend(bg="white",bty = "n","bottomright",paste("n = ",floor(tt$parameter/1000),",", 2+tt$parameter%%1000,sep=""),cex=1.5)
lines(lowess(x=agg67d2$dif_log_mean_amount,y=agg67d2$dif_log_mean_hours_cr0,f=0.68),lty=2)
#Panel 2
tt<-cor.test(agg6d2$dif_log_mean_amount,agg6d2$dif_log_mean_hours_cr0)
plot(main=c("(2) Contract Lobbyists"),cex.main=1.5,cex.lab=1.2,cex.axis=1.5,x=agg6d2$dif_log_mean_amount,y=agg6d2$dif_log_mean_hours_cr0,
     xlim=c(-8,2),ylim=c(-8,2),cex=1,col="gray",xlab="Ln(Dyad Amount)-Ln(Dyad Biennium Mean Amount)",ylab="Ln(Dyad Hours)-Ln(Dyad Biennium Mean Hours)")
text(x=-5.75,y=1.85,paste("r = ",round(tt$estimate,3)," (",round(sqrt( (1-tt$estimate^2)/(tt$parameter)),3),")",sep=""),cex=1.5)
legend(bg="white",bty = "n","bottomright",paste("n = ",floor(tt$parameter/1000),",", 2+tt$parameter%%1000,sep=""),cex=1.5)
lines(lowess(x=agg6d2$dif_log_mean_amount,y=agg6d2$dif_log_mean_hours_cr0,f=0.68) ,lty=2)
#Panel 3
tt<-cor.test(agg7d2$dif_log_mean_amount,agg7d2$dif_log_mean_hours_cr0)
plot(main=c("(3) In-House Lobbyists"),cex.main=1.5,cex.lab=1.2,cex.axis=1.5,x=agg7d2$dif_log_mean_amount,y=agg7d2$dif_log_mean_hours_cr0,
     xlim=c(-8,2),ylim=c(-8,2),cex=1,col="gray",xlab="Ln(Dyad Amount)-Ln(Dyad Biennium Mean Amount)",ylab="Ln(Dyad Hours)-Ln(Dyad Biennium Mean Hours)")
text(x=-5.75,y=1.85,paste("r = ",round(tt$estimate,3)," (",round(sqrt( (1-tt$estimate^2)/(tt$parameter)),3),")",sep=""),cex=1.5)
legend(bg="white",bty = "n","bottomright",paste("n = ",floor(tt$parameter/1000),",", 2+tt$parameter%%1000,sep=""),cex=1.5)
lines(lowess(x=agg7d2$dif_log_mean_amount,y=agg7d2$dif_log_mean_hours_cr0,f=0.68) ,lty=2)
#Panel 4
tt<-cor.test(agg67d3$dif_log_mean_amount,agg67d3$dif_log_mean_hours_cr0)
plot(main=c("(4) All Lobbyists - Yearly Aggregation"),cex.main=1.5,cex.lab=1.2,cex.axis=1.5,x=agg67d3$dif_log_mean_amount,y=agg67d3$dif_log_mean_hours_cr0,
     xlim=c(-8,2),ylim=c(-8,2),cex=1,col="gray",xlab="Ln(Dyad Amount)-Ln(Dyad Mean Amount)",ylab="Ln(Dyad Hours)-Ln(Dyad Mean Hours)")
text(x=-5.75,y=1.85,paste("r = ",round(tt$estimate,3)," (",round(sqrt( (1-tt$estimate^2)/(tt$parameter)),3),")",sep=""),cex=1.5)
legend(bg="white",bty = "n","bottomright",paste("n = ",floor(tt$parameter/1000),",", 2+tt$parameter%%1000,sep=""),cex=1.5)
lines(lowess(x=agg67d3$dif_log_mean_amount,y=agg67d3$dif_log_mean_hours_cr0,f=0.68) ,lty=2)
#Panel 5
tt<-cor.test(agg6d3$dif_log_mean_amount,agg6d3$dif_log_mean_hours_cr0)
plot(main=c("(5) Contract Lobbyists - Yearly Aggregation"),cex.main=1.5,cex.lab=1.2,cex.axis=1.5,x=agg6d3$dif_log_mean_amount,y=agg6d3$dif_log_mean_hours_cr0,
     xlim=c(-8,2),ylim=c(-8,2),cex=1,col="gray",xlab="Ln(Dyad Amount)-Ln(Dyad Mean Amount)",ylab="Ln(Dyad Hours)-Ln(Dyad Mean Hours)")
text(x=-5.75,y=1.85,paste("r = ",round(tt$estimate,3)," (",round(sqrt( (1-tt$estimate^2)/(tt$parameter)),3),")",sep=""),cex=1.5)
legend(bg="white",bty = "n","bottomright",paste("n = ",floor(tt$parameter/1000),",", 2+tt$parameter%%1000,sep=""),cex=1.5)
lines(lowess(x=agg6d3$dif_log_mean_amount,y=agg6d3$dif_log_mean_hours_cr0,f=0.68) ,lty=2)
#Panel 6
tt<-cor.test(agg7d3$dif_log_mean_amount,agg7d3$dif_log_mean_hours_cr0)
plot(main=c("(6) In-House Lobbyists - Yearly Aggregation"),cex.main=1.5,cex.lab=1.2,cex.axis=1.5,x=agg7d3$dif_log_mean_amount,y=agg7d3$dif_log_mean_hours_cr0,
     xlim=c(-8,2),ylim=c(-8,2),cex=1,col="gray",xlab="Ln(Dyad Amount)-Ln(Dyad Mean Amount)",ylab="Ln(Dyad Hours)-Ln(Dyad Mean Hours)")
text(x=-5.75,y=1.85,paste("r = ",round(tt$estimate,3)," (",round(sqrt( (1-tt$estimate^2)/(tt$parameter)),3),")",sep=""),cex=1.5)
legend(bg="white",bty = "n","bottomright",paste("n = ",floor(tt$parameter/1000),",", 2+tt$parameter%%1000,sep=""),cex=1.5)
lines(lowess(x=agg7d3$dif_log_mean_amount,y=agg7d3$dif_log_mean_hours_cr0,f=0.68) ,lty=2)

